Assessing the impact of land use land cover change on regulatory ecosystem services of subtropical scrub forest, Soan Valley Pakistan

This study investigated the effect of land use land cover (LULC) changes on carbon sequestration in the Hayat-ul-Mir subtropical scrub reserve forest, Pakistan. Supervised maximum likelihood classification of Landsat satellite imagery was done to assess spatio-temporal changes in LULC during 2007, 2013 and 2019. The CA–Markov model was used to simulate LULC of 2030. Spatial LULC data and carbon pools data was processed in Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) carbon model to investigate the effect of LULC on future carbon dynamics. The analysis revealed increase in cover of A. modesta and O. ferruginea and decrease in agriculture, built up and barren area of forest during 2007–2019 and 2030. The analysis also showed that the forest would additionally sequester 111 Mg C with an overall Net Present Value of $4112.05 in year 2030. The analysis revealed LULC changes on 25% area with increase and decrease in the value of ecosystem service (at some location) from carbon storage and loss as CO2 emissions respectively depending on the type of LULC converted. The study is helpful in identifying areas of potential carbon sequestration to maximize net benefits from management interventions.

NDVI classification. NDVI (Normalized difference vegetation index) time series maps (year 2007, 2013 and 2019) were developed by measuring the reflectance of red and infrared bands of satellite imagery based on density and intensity of vegetation using Eq. (1) 23 . The imagery was downloaded from United States Geographical Survey (USGS-Earth Explorer) website and was processed in ArcMap (v. 10.2). The reference years were selected to analyze the state of forest in the pre and post implementation phase of the Forest Act (Punjab Forest [Amendment] Act 2010 in case of this study). To compare NDVI values, all satellite images were acquired for the same month of the respective years (i.e. October) after monsoon rains when vegetation was stable and also to avoid seasonal variations ( Table 1). The difference in the dates of image acquisition (although in the same month) was unavoidable because of cloud cover in the satellite imagery and image distortion. It was also based on availability of satellite data for a particular date in the month of October.
where NIR and RED are the reflectance in near infrared and red bands respectively based on the density of vegetation. The NDVI values generally range between − 1.0 and + 1.0 referring to no or high (healthy) vegetation cover. A value of 0.1 or below indicates bare soil or rock, however, value of 0.2 to 0.3 represents shrub/ low vegetation cover and higher values represent trees/high vegetation cover 24 .
Valuation of ecosystem carbon storage and sequestration services. Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) carbon model (v. 3.3.3) was used for the quantification, mapping and valuation of the effect of LULC changes on ecosystem carbon storage and sequestration. The model was fed with LULC spatial data for year 2019 and carbon data from field estimates of four pools i.e. aboveground, belowground, soil and litter (dead) to estimate carbon stored in each grid cell. Carbon stock of the entire ecosystem www.nature.com/scientificreports/ was estimated by aggregating individual carbon estimates of each LULC class (Eqs. (2), (3)). Future LULC raster (for year 2030) was input into the model to predict carbon sequestration (Eq. (4)) in future 25 . For economic valuation, the country specific median price of a $36.67/ton of carbon (the social cost of carbon) with a discount rate of 2% was used 26 . The model resulted in three dimensional estimate of carbon stored and sequestered in each LULC class along with its monetary benefit.
where A is the area of grid cell (0.09 ha), Ca m,i,j , Cb m,i,j , Cs m,i,j and Cd m,i,j refers to the aboveground, belowground, soil and dead (litter) carbon density (Mg ha −1 ) for grid cell (i, j) for land use type (m). C is the aggregate carbon density of the ecosystem for the current year (C cur ) and predicted year (C fut ) and S is the carbon sequestration across the study area 25 . InVEST input data preparation. Land Cm, i, j,   (6)). Cellular Automata (CA) model was used for spatial simulations. In CA, each cell decided the data between different states in time for itself and its neighboring cells (through contiguity filter) and decided change in each cell by rules. The weight factor is applied based on proximity of nuclear and neighbor cells which is then combined with transition probabilities to predict the state of closer cells so that land use change is not a random decision. A combined CA-Markov model was used to simulate the spatio-temporal dynamics among the LULC classes using transition probability matrix as input to CA model to predict future LULC 11,25 . where S is the status of land use at time t and t + 1; and Pij is the transition probability matrix.
Using CA-Markov model in the IDRISI Selva (v. 17.0), number of CA iterations was set to 11 as time interval, 5 × 5 contiguity filter and transition probability matrix (2013 to 2019) was used to predict LULC of 2030. Before running the final model, it was validated by creating a prediction map of 2019 from 2007 and 2013 LULC data and then it was compared with the actual (real) LULC map of 2019 using kappa statistics (Eq. (7)). The LULC maps generated were visualized in ArcGIS (v. 10.2), study area was extracted and converted to InVEST compatible file format before running sequestration and valuation model. The LULC change map was developed in ArcMap (v. 10.2) to visualize inter-conversion of LULC classes during 2019 to 2030. Complete scheme of study is presented in Fig. 2.
where kappa is index to test simulation accuracy; Po, Pc and Pp are the actual, expected and ideal (100%) simulation accuracies respectively 11 .
Assessment of forest carbon pools. Estimates of the above and belowground carbon stored in A. modesta and O. ferruginea were adopted from our recent study 19 conducted in the HM forest by measuring tree dendrometric variables in 47 survey plots (0.04 ha). Sub-plots of 1 m 2 were marked to collect leaf litter (dead fraction) and soil samples (at 15 cm depth). To determine carbon stored in the dead pool, 1 g of oven dried leaf litter was ignited at 550 °C in a muffle furnace for 2 h following Salehi et al. 27 . Soil total organic carbon (TOC) was determined through Walkley Black chromic acid wet oxidation method after reducing organic compounds of soil sub-sample (1 g) with 1 N K 2 Cr 2 O 7 . The unreduced dichromate was determined by oxidation-reduction titration. TOC (%) was multiplied with sample depth and soil bulk density to estimate soil organic carbon stock (SOC) on per hectare basis 28 .

Results
NDVI cover. Figure   www.nature.com/scientificreports/ ground carbon (AGC), belowground carbon (BGC), soil and dead carbon. The limited agriculture activity on the forest boundary was recorded with storing 3.8 ± 0.9 Mg/ha of carbon whereas least carbon was stored as SOC in the built-up and barren LULC classes (Table 4). Multiple linear regression analysis revealed a significant positive relationship of LULC with total carbon stored (p < 0.05) in the pools of each class i.e. AGC, GBC, SOC and dead C. The analysis suggests increase in carbon storage with increase in the vegetation density of LULC class. Similarly, positive significant relationship was observed between LULC class and corresponding area covered (p < 0.001) ( Table 5).
The InVEST carbon model helped in identifying the LULC class with higher carbon storage potential. The tool estimated maximum 15.22 Mg/ha of total carbon stored in the vegetative cover of the forest in the base (2019) and predicted (2030) land covers. However, least carbon was stored (0.33 Mg/ha) in agriculture, built-up and barren land cover classes suggesting their least contribution in the provision of ecosystem service (Fig. 5a,b). Changes or interconversions in LULC classes between base (2019) and predicted (2030) year were detected to determine their benefit and potential for carbon sequestration (Fig. 5c-e)

Discussion
Current study presents modeling based prediction of the expected LULC changes on regulatory ecosystem service of the HM forest on a spatio-temporal scale. The study site is under jurisdiction of provincial government and thus is a declared reserve forest under Punjab Forest (Amendment) Act 2010 prohibiting any cultivation, www.nature.com/scientificreports/ land clearings, grazing, mining and timber harvesting unless granted under limited rights 29 . Therefore, limited human activity was observed in the HM reserve forest except some agricultural activity and built up area on forest edges only. Barren land or area with no vegetative cover was seen throughout the forest, but it declined and was replaced with tree cover overtime (2007 to 2019) as also represented through NDVI.  www.nature.com/scientificreports/ NDVI helped in investigating density and health of the forest on spatio-temporal basis. The most dominant cover class in the HM forest was made by the bi-climax community of O. ferruginea and A. modesta represented with an increase in NDVI value overtime. The limited human activity in HM forest reserve resulted in increase (~ 32%) in its vegetation density and hence biomass carbon stock during 2007 to 2019. Reason for this increase in the vegetation density can be attributed to the implementation of Punjab Forest (Amendment) Act 2010. The act defined rights and restricted human interference in the form of encroachments, construction, timber and fuelwood harvesting, damage to soil and natural vegetation, cattle grazing, trespassing and any kind of land use change other than for the conservation purpose under sub-section 1 of section 26. The act also revised penalties to offences in the form of imprisonment for six months or fine ranging between PKR 10,000 to PKR 0.5 million or both depending on the worth of damage which could double if offence committed after sunset, as defined under sub-section 2-4 of section 26 29 . Consistent to the findings of present study, Kamwi et al. 30 also reported increase in the forest cover of protected and communal areas of Zambezi region in Namibia during 1991-2010 compared to 1984-1991 after the implementation of forest policy restricting any human interference. In the absence of forest policy and participatory forest management, 37.38% decline in the forest area was observed in 2019 compared to 1991 (20%) due to agricultural expansions (+ 33.5%) in Komto protected forest of Ethiopia which is contrary to the findings of current study 31 .
Other reason for increase in the vegetated cover of HM forest is related to a significant reduction in the livestock grazing during the last 10 years due to shift to non-livestock dependent livelihood practices of the adjacent communities. This is evident from a significant increase in the number of factories from just 28 in year 2007 to 145 in year 2019 in the district. There was also a significant increase in the available agriculture machinery which reduced dependence on work animals [32][33][34] . The decrease in grazing pressure provided relief to the scrub forest vegetation (both shrubs and climax species) and it restored as reflected in the NDVI and LULC for year 2013 and 2019. Findings of this study are contrary to Mannan et al. 21 reporting decrease in vegetation cover with lower NDVI values (1998 to 2018) for subtropical and moist temperate forests dominated with A. modesta, O. ferruginea, Ziziphus spp. and Pinus spp. on the foothills of the Himalayas, Pakistan due to increase in urbanization and human interferences. Zhao et al. 11  Apart from the LULC changes, other factors such as rising atmospheric CO 2 concentration at a rate of 2.5 ppm/year (expected 438 ppm in 2030) could induce a fertilization effect enhancing the photosynthetic rate and supporting the growth of both shrub and climax species 35 . But, coupled with the elevated CO 2 levels, future warming and change in precipitation could limit the stem and root growth with increase in soil respiration 36 . Ghafoor et al. 37 reported significant reduction in the regeneration rate and growth of A. modesta and O. ferruginea in HM subtropical forest under short term climate warming. Considering the current pace of global warming, where already 1 °C rise in temperature has been observed over the last 40 years in Pakistan, decrease in precipitation in the Soan valley and future warmer climate with elevated CO 2 concentration might result in decrease in ecosystem productivity 38,39 . But the model was limited in simulating the effect of climatic variability on land cover of the forest.
Gradual increase in tree cover during 2007-2019 and simulation for 2030 implies accumulation of biomass and higher carbon sequestration potential of the HM forest in the future. The order of carbon stored in the HM forest was as; aboveground > soil > belowground > dead pools. The aboveground vegetation was identified as the largest pool storing carbon where major contribution was made by O. ferruginea, but comparatively the estimates were lower than other subtropical forests of Pakistan as reported by Nizami (2012) and Shaheen et al. 40 . The  www.nature.com/scientificreports/ estimates of soil carbon density of the current study are higher than those reported by Nizami 17 and represent its significance in storing carbon for long term if undisturbed. The biophysical assessment of tree carbon stock showed that the subtropical broadleaved evergreen species (O. ferruginea) stored more carbon than subtropical broadleaved deciduous species (A. modesta) in the HM forest reserve. This interspecific difference suggests the relative importance and potential of both species in carbon sequestration as indicated through the InVEST carbon model. A similar trend in relative carbon storage and sequestration was reported by Arunyawat and Shrestha 41 for deciduous and evergreen plantations in Northern Thailand. The ecosystem service assessment presented an increase in carbon storage in the HM forest in the future under BAU scenario which is consistent to the findings of Sil et al. 8 and Zhao et al. 11 for the simulated period in the river basins in Portugal and China respectively. Carbon storage was found comparatively high at the sites of evergreen species than where deciduous species is distributed in the HM forest, whereas least carbon storage in the non-vegetated classes in the base and predicted year. The historical trend in the relative area coverage of both species during 2007-2019 and predicted LULC presented O. ferruginea dominating the A. modesta suggesting that under a future warmer climate the evergreen species will have more potential in mitigating climate change 37 . This analysis is also helpful in identifying areas of potential carbon storage and where to implement management interventions to increase the service of tree carbon sequestration. www.nature.com/scientificreports/ The land use change analysis reported a large area of the forest (75%) where no change in LULC is expected during 2019 to 2030. The simplified carbon cycle in the InVEST carbon model assumed a linear change in carbon sequestration overtime and did not calculate SR where no change in LULC is expected to occur during base and predicted year. Therefore, the model resulted in a low sequestration estimate of 111.91 Mg C (or 0.006 Mg C/ ha/year) during 2019 to 2030 for the remaining land area (25%) which is expected to undergo LULC change. Consistent to the present study, Zhao et al. 11 also reported low sequestration rates due to model limitations.
The analysis also identified sites with negative sequestration rates where loss of stored carbon (as CO 2 emissions) is expected to occur due to conversion of vegetated land to built up or barren area in future. Additionally, low SR was calculated for the land cover class where area under O. ferruginea is expected to be replaced by A. modesta on lower elevations. The low SR for this conversion is related to a relatively low potential of A. modesta (deciduous) to store carbon than O. ferruginea (evergreen) and InVEST assumed it as some loss of stored carbon. Increase in ecosystem service was observed where A. modesta is expected to be replaced by O. ferruginea on higher elevations in the future and it was highest for the pixels where agriculture, built up and barren area are expected to convert to forest vegetation.
Considering the intensity of LULC change in the HM forest, SR was low compared to the other studies. Sil et al. 8 reported higher SR (1.63 Mg/ha/year) under projected (2020) forest expansion scenario in the Sabor river basin, Portugal due to 50% increase in broadleaved, chestnut and coniferous forests. The difference in the estimates of both studies is related to the type and extent of LULC converted and site climatic factors affecting plant growth. Compared to this study, Zhao et al. 11 predicted 16% of the shrub land of the Heihe river basin in China to be converted into grassland and bare land in 2029 reducing the net carbon density. Leh et al. 12 also documented decrease in ecosystem carbon storage and sequestration service during 2000-2009 from LULC changes in West Africa. Carbon density in a forest does not remain constant, but it fluctuates over the life cycle of a species with high sequestration ability during its full vigor growth phase and falls with tree age and thus also affects carbon additions to soil from litter production. Increase in site productivity by implementation of management interventions also improves ecosystem potential of carbon storage and sequestration 42 . Thus the oversimplification of the carbon cycle in InVEST model and assumption of zero SR where no LULC change is observed might lead to some misrepresentation of the actual supply of this ecosystem service.
Like carbon SR, gains and losses in the economic value of ecosystem service was observed for the type of LULC converted during 2019-2030. The monetary benefit of carbon sequestration is expected to be highest at the sites where non-vegetated cover classes might convert into O. ferruginea and A. modesta and decrease in it where the forest cover is expected to be replaced with barren land. The net present value (NPV) of carbon sequestered in the current and future land covers ($4112 or PKR. 0.66 million only) of the HM forest was low because of LULC change in 25% of the land area only and limitation of the model to calculate NPV on remaining 75% of the intact area. The monetary value of this ecosystem service is difficult to be compared with other studies due to the difference in currency units, the social cost of carbon (SCC) and discount rate reflecting society's preference for immediate benefits of carbon sequestration over future benefits. Contrary to the findings of this study, Sil et al. 8 with SCC of $ 23, $ 44 and $ 312 (at a discount rate of 1, 3 and 7%) reported a significantly high value of this ecosystem service from approximately 50% change in LULC under future scenario. Together with the estimates of carbon storage and sequestration, the economic valuation of this service can help in identifying the areas of potential carbon sequestration to maximize the net benefits of forest management interventions.

Conclusion
The current study reveals that land cover changes have occurred in the Hayat-ul-Mir forest during past decade. Over time (2007 to 2019), the land covers of A. modesta and O. ferruginea were increased and are further expected to increase in the simulated period (year 2030) based on BAU scenario. The current forest management practices and the restrictions imposed by the act were found satisfactory in reducing human impact although limited agricultural activity on the forest edges was observed but in the declining trend. Built up and the barren land area covers also decreased overtime and most of these classes were represented around coal mining sites and shrine during current and simulated land covers. The biophysical and economic assessment of the service suggested increase in future carbon storage with current management practices. The InVEST model was limited in the assessment of carbon sequestration and did not calculate it for the area which is expected to remain same in 2030. This might had led to a misrepresentation of actual carbon sequestration potential and net benefits of the HM forest. Improvements in the model are therefore required to accurately account for the carbon sequestration rates. Some interconversions of the LULC classes were also observed for the simulated land cover where the vegetation might replace the other land uses. For those parcels of land, such interconversions would bring positive impact in terms of carbon storage. The model estimated negative rates for areas where termination of agriculture and mining activity (representing some built up area) would make land vacant. CO 2 emissions are expected from those land parcels if not restored with native vegetation. Hence based on this analysis, the study provides a valuable input to inspect the land uses that are worse off and where appropriate management interventions must be introduced to increase the carbon storage and sequestration potential of the subtropical forest.

Data availability
The data generated and analyzed in this study is included in this article.